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1.  Introduction 


The  collocation  method  is  now  widely  accepted  and  has  been  used  to 
solve  many  geodetic  problems,  e.  g.  the  computation  of  geoidal  undulations  and 
deflections  of  the  vertical  (Tscherning,  1973;  Tscherning  and  Rapp,  1974; 
Lachapelle  1975a  and  1975b)  as  well  as  for  the  extension  of  the  gravity  field 
from  satellite  altimetry  (Rapp,  1974  and  1976  ).  In  its  narrower  and  original 
form  this  technique  has  been  used  for  the  prediction  of  point  or  mean  anomalies 
from  the  known  gravity  field  (Rapp,  1964  ).  hi  fact,  the  method  has  been  used, 
at  least  experimentally,  to  assist  in  the  solution  of  all  aspects  of  the  geodetic 
problem. 


The  covariance  function  is  basic  to  collocation  methods.  This  fact 
is  emphasized  by  Moritz  in  his  important  work  on  the  theory  of  collocation  when 
he  comments  that,  "Not  only  conceptually,  but  also  computationally  the  covariances 
represent  the  crucial  point  in  the  present  method. . . " (Moritz,  1972,  p.  21). 

These  functions  are  meant  to  represent  the  statistical  characteristics  of  the 
fields  they  represent,  and  in  their  derivation  conditions  of  homogeneity  and 
isotropy  are  assumed  to  exist  in  the  fields.  It  is  the  correctness  of  this 
assumption  which  is  being  investigated  here;  if  the  assumptions  are  found  to 
be  inaccurate  it  will  be  necessary  to  assess  the  impact  that  the  lack  of  accuracy 
has  on  the  method,  and  to  suggest  ways  of  modifying  the  method  to  account  for 
an  'improved'  model. 

Although  the  covariance  functions  of  the  anomalous  potential  field  is  the 
theoretical  basis  for  covariance  functions  of  other  (potential-related)  fields 
(Jordan,  1972;  Moritz,  1972,  p.  94-99)  it  is  usual  for  the  Ag  field  to  form  the 
basis  for  most  collocation  solutions  (e.g.  Lachapelle,  1975a,  p.  23-24).  For 
this  reason  this  investigation  will  at  first  concentrate  on  the  detection  of  anis- 
otropy in  anomalous  gravity  fields.  This  is  then  extended  to  the  behavior  of 
cross -covariance  functions  and  the  effects  of  using  an  isotropic  model  for 
collocation  in  an  anisotropic  region  are  investigated. 

The  assumption  of  an  homogeneous  gravity  field  throughout  the  globe 
has  been  queried  for  as  long  as  covariance  techniques  have  been  discussed.  For 
example,  Rapp  investigated  three  different  regions  of  the  United  States  and  found 
significantly  different  covariance  functions  for  each  from  the  analysis  of  these 
fields  (Rapp,  1964,  p.42-65).  Williamson  and  Gaposbkin  analyzed  subsets  of 
data  of  ocean  and  continental  gravity  assuming  isotropy,  and  found  the  differences 
between  the  covariances  to  be  significant.  They  concluded  that  ". . .gravity  is 
not  stationary"  (Gapoahkin,  1973,  p.208).  Similarly,  in  a statistical  analysis  of 
anomalous  gravity  of  Czechoslovakia,  Vyskoiil  (1970)  investigated  the  covariance 
functions  for  Ag  means  of  10'  x 15'  blocks  and  concluded  "the  theoretical 
assumption  of  homogeneity  and  Isotropy. . . is  basically  only  forced  and  in  reality 
it  will  probably  not  always  be  satisfied.  " (See  also  Kaula,  1966b,  p.  58.) 
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It  is  common  practice  to  allow  for  the  non-homogeneity  of  the  gravity 
field  by  restricting  the  analysis  to  that  data  pertinent  to  the  area  of  solution. 

That  is,  a local  covariance  function  based  on  a restricted  gravity  field  is  derived 
and  used  in  subsequent  computations  for  that  area.  The  problem  remaining  is 
how  to  test  such  a local  field  for  isotropy  and,  if  characteristics  of  non-isotropy 
are  exhibited,  how  to  include  these  in  the  computations. 

As  already  mentioned,  Vyskofiil  (qv)  has  expressed  concern  that  non- 
isotropy  can  exist.  A comparison  of  the  profile  analysis  by  Rapp  (1964,  p.44-56) 
shows  large  differences  between  N-S  and  E-W  profiles  of  the  same  area,  casting 
doubt  on  the  assumption  that  the  covariance  functions  are  independent  of  azimuth. 
Kubifckovi  (1974)  has  studied  the  gravity  field  in  the  Carpathian  Mountains  of 
Czechoslovakia,  and  concludes  that  the  departure  from  a situation  of  isotropy  is 
not  insignificant,  suggesting  that  use  of  a covariance  function  based  on  an  isotropic 
model  "is  impermissible  without  prior  analysis  of  the  consequences. " 

The  purpose  of  this  report  is  threefold.  Firstly,  it  is  to  see  how 
anisotropic  characteristics  in  a data  field  are  best  detected  and  portrayed. 
Secondly,  to  find  out  if  an  anisotropic  model  can  be  introduced  into  the  theory  of 
prediction  and  collocation.  Lastly,  we  wish  to  investigate  if  an  'improved' 
statistical  model  produces  any  significant  effect  on  the  results  of  prediction  and 
collocation  solutions. 

It  may  be  that  such  improvement  is  marginal.  It  has  already  been 
pointed  out  that  ". . .the  choice  of  the  covariance  function  is  less  critical  (i. e. 
than  internal  consistency  in  the  computations  - author's  note)  because  the  results 
of  least  squares  collocation  are  not  very  sensitive  with  respect  to  the  covariance 
function  chosen,  in  the  same  way  that  the  results  of  ordinary  least  squares 
adjustment  do  not  depend  strongly  on  the  weights"  (Moritz,  1972,  p.  90).  Tukey 
(1970,  p.  165-166)  estimates  that  in  ordinary  least  squares  situations,  if  the 
ratio  of  each  true  weight  to  the  weight  you  use  does  not  vary  by  more  than  a 
factor  of  2,  then  the  efficiency  of  the  fit  is  always  at  least  88. 8%.  However,  he 
then  differentiates  between  the  normal  least  squares  solution  and  the  time-series 
situation,  pointing  out  that  due  to  the  overlap  in  data,  some  bias  will  result  if  the 
variances  and  covariances  are  incorrectly  estimated. 

"When  in  doubt,  smooth"  (Moritz,  1972,  p.  125)  is  a valuable  and 
useful  maxim,  but  there  may  be  cases  where  such  a procedure  is  contrary  to  a 
preliminary  appreciation  of  the  situation.  It  may  be  that  to  'smooth'  without 
fully  establishing  the  'doubt'  is  detrimental  to  the  resultant  solution. 
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2.0  The  Theoretical  Background 

2. 1 The  Covariance  Function 


The  application  of  covariance  analysis  to  gravimetric  geodesy  is  well 
known  (Heiskanen  and  Moritz,  1967,  252-259,  266-274).  The  fundamentals  will 
be  briefly  reviewed  here  in  order  to  emphasize  how  the  assumption  of  random- 
ness of  the  gravity  field  influences  the  mathematical  representation  of  Its 
covariance  function. 


'The  covariance  characterizes  the  statistical  correlation  of  the  two 
signals,  Ag  and  Ag* , which  is  their  tendency  to  have  the  same  size  and  sign. 

If  the  covariance  is  zero  then  the  anomalies  Ag  and  Ag*  are  uncorrelated 

i.e.  the  size  and  sign  of  Ag  has  no  Influence  on  the  size  and  sign  of  Ag*"(Ibid,  253). 


This  correlation  is  estimated  by  taking  the  mean  of  all  product  pairs 
throughout  the  globe,  thus: 


i r2TT  rn  r2TT 

(2.1)  C(P,Q)  = —p  J^  = 0 J0  = o J a=0  Ag  (^,A»)  Ag(^,Xq)  sin0  d8  dX  da 


where 

P,  Q represent  the  current  points  in  the  product  pair,  apart 
0,  X are  the  spherical  coordinates  of  these  points,  and 
a is  the  azimuth  from  P to  Q. 


If  isotropy  is  assumed  the  covariance  is  independent  of  azimuth,  and 

(2.1)  becomes: 


(2.2) 


£ Lo  Lo  ** su>  9 d9<u 


If  homogeneity  is  assumed  the  covariance  is  independent  of  location  on 
the  globe,  and  since: 


(2.3) 


.2tt  „tt 


£ Lo  Je-0  £//«• 


(2. 2)  becomes  a function  of  separation  only.  Thus: 


(2.4) 


C<P,  Q)  - J*J  Agq  da 
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If  local  regions  are  considered,  the  reference  model  can  be  assumed 
planar  and  the  separation  expressed  as  the  straight  line  distance  r*j , where: 


(2.  5) 


C(r)  =*  M 


{Agp  Agqj- 


Cpq 


M | x | being  the  mean  value  of  x. 

For  r = 0 

(2.6)  Co=C(0)  = M {Ag2} 


the  variance  of  the  gravity  anomalies.  Note  also  that: 
(2.7)  . M {Ag}  = 0 


which,  if  not  the  case,  must  be  enforced  by  "centiring"  the  data  before  commencing 
the  analysis  (see  Rapp,  1964,  p.7-8). 

In  a local  area  the  means  of  product  pairs  are  usually  computed  from 
point  data  or  from  means  of  small  areas  (e.  g.  0°1  x 0®  1)  arranged  on  a grid  pattern. 
If  the  field  is  isotropic,  the  values  of  the  means  for  any  one  separation  are 
independent  of  azimuth.  Conversely,  if  the  means  are  computed  without  reference 
to  azimuth  (see  equation  2.2)  the  property  of  isotropy  is  enforced  onto  the  field. 

An  analytical  function  can  be  fitted  to  the  values  of  C (r)  (e.  g.  Moritz, 
1976,  p.  25  - 29).  This  ftmction  is  a convenient  way  to  generate  the  elements  of  the 
matrices  used  in  subsequent  computations.  However,  this  fitting  will  have  probably 
introduced  further  "smoothing"  into  the  covariance  function;  thus  the  price  paid  for 
convenience  may  be  a loss  of  precision  in  the  representation  of  the  statistical 
relationships  between  two  points  in  the  subject  data  field. 

A review  of  the  process  will  underline  the  stages  at  which  the  actual 
covariance  may  have  undergone  change.  The  covariance  function  has  been  found 
by: 

(I)  taking  the  means  of  product  pairs  throughout  the  region  being 
analyzed, 

(II)  taking  the  means  regardless  of  the  directional  relationships  of 
successive  point  pairs,  and 
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(ill)  fitting  an  analytical  function  to  the  resultant  discrete  field. 


If  homogeneity  is  present,  the  smoothing  introduced  by  (i)  will  be 
minimal.  If  isotropy  exists,  step  (ii)  will  also  introduce  a minimum  of  smoothing. 
If  both  conditions  are  present,  then  the  function  obtained  in  (iii)  will  accurately 
reflect  the  existing  situation.  How  the  covariance  function  is  used  will  be 
demonstrated  in  the  following  section. 


2.2  Least  Squares  Prediction 


This  discussion  will  concentrate  on  the  application  of  stochastic  methods 
in  the  prediction  of  gravity  anomalies.  This  will  assist  in  the  appreciation  of  the 
stochastic  techniques  generally,  and  in  particular  will  underline  how  the  assumption 
of  isotropy  in  die  data  field  carries  through  into  the  formulation  and  solution  of  the 
problem.  The  more  involved  multivariate  case  will  be  treated  separately  (see 
section  2.4). 

The  covariance  relationship  is  used  in  the  determination  of  the 
coefficients  apl  of  a univariate  Kolmogorov-Wiener  prediction  model  (Grafarend 
and  Offermans,  1975,  3;  Grafarend,  1976,  p.152).  Thus: 

(2.8)  Ag'p  = Y Ag,  1 = 1,  n 

— J 

l 

Ag > being  the  predicted  value  of  Ag  at  P. 

The  error  of  prediction  is  found,  thus: 


(2.9) 


These 

(2. 10) 


O = ( Agp  - V ap,  Ag  ,)a 
1 

squared  errors  are  now  averaged  through  the  region  for  all  points  P.  Thus: 

M = M | Agp3 j-  - 2V  ap,  M -[Agp  Agij-  + V V a*  ap,  M-j^Ag,  Agjj 

1 1 J 

l ■ l,n 
J»l,n 
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Assuming  that: 


(2.11a)  M-[Agfaj-  = 

(2.11b)  M-j^Agr  Agtj-  = M-^AgtAgj ]■  for  ry=  r>t 

and  using  the  covariance  function  derived  in  section  (2. 1),  (2. 10)  becomes: 
(2.12)  njf3»  Co  - 2Y  art  Cn  + YY  an  aM  Cm 

— I -J  _» 

1 1 1 


where 

(2.13)  m?  = M aj- 


Note  the  approximations  in  (2.11).  The  values  for  Agr  are  (obviously) 
unknown,  and  so  some  assumption  regarding  their  statistical  relationship  with 
respect  to  known  values  ( Agj)  must  be  made  for  (2. 10)  to  be  solveable.  The 
accuracy  of  the  estimates  made  in  (2. 11)  will  depend  on: 


(i)  how  well  the  characteristics  of  homogeneity  and  isotropy  describe 
the  actual  Ag  field,  and 

(ii)  the  amount  of  data  used  to  find  the  C0  , Cn  and  values. 


This  second  point  is  fairly  obvious  and  not  of  direct  concern  at  the 
moment.  It  is  sufficient  to  say  that  if  the  data  is  sparse,  then  a poor  estimate 
for  the  covariance  function  is  likely.  This  is  despite  any  indications  of  good 
accuracy  obtained  from  an  estimate  of  internal  precision  (see  section  2. 3). 

The  first  point  is  of  more  importance  in  this  present  discussion.  If 
the  field  is  well  modelled  by  a stationary  statistical  model,  then  the  approxima- 
tions for  M^Ag  a|  and  M ^Agr  AgtJ  will  be  realistic.  Otherwise,  these 

approximations  are  "smoothing"  the  real  situation  and  characteristics  of  homo- 
geneity and  Isotropy  are  Imposed  onto  the  predicted  field. 


The  coefficients  a»i  in  equation  (8)  are  now  solved  in  the  usual  way  by 
least  squares,  viz. 


o. 

d a*i 


Thus  it  is  found  that: 


(2.14)  a,;  » £ Cj1  C,}  J - 1,  n . 

i 

The  predicted  value  is  found  by  substitution  of  (2. 14)  into  (2. 8).  Thus 
it  can  be  seen  how  the  adopted  statistical  model  is  imposed  (by  way  of  Cu  and 
Cm)  onto  the  predicted  values  for  Agr  • 


2. 3 Accuracy  of  Least  Squares  Prediction 


The  investigation  is  carried  into  the  estimation  of  accuracy  to  help 
emphasize  the  role  which  the  statistical  model  plays  in  this  phase  of  the  develop- 
ment. 


Equation  (2.9)  can  be  used  as  the  starting  point;  thus: 


(2.15) 


- y 


an 


where  a?t  are  found  from  (2. 14).  In  equation  (2. 15)  we  are  comparing  the  pre- 
dicted value  of  Agr  with  Agr  assumed  known. 

The  procedure  described  in  section  (2. 2)  to  find  the  mean  square  errors 
is  again  followed  and  it  is  found  that: 


(2. 16) 


m * - C0  - Cn  Cjl  . 


(In  (2. 16)  matrix  notation  is  used  to  simplify  the  expression).  It  must  be 
noticed  that  estimates  (2. 11a)  and  (2.  lib)  are  again  used.  The  errors  being 
determined  are  therefore  not  those  at  the  unknown  points  P but  at  known  points  i. 
That  is,  equation  (2. 16)  is  a measure  of  the  "fit”  of  a stationary  statistical  model 
onto  the  known  data.  (See  also  Grafitrend,  1976,  p.  154.) 

This  is  an  Important  distinction  to  make.  While  m * is  probably  a 
useful  guide  to  the  accuracy  of  the  predicted  values,  it  is  strictly  a measure  of 


the  "fitness  of  the  model".  As  such,  it  will  be  useful  when  comparing  predictions 
based  on  anisotropic  models  with  those  using  isotropic  models  (see  section  4. 1). 

From  the  above  comments  one  can  see  the  danger  in  placing  too  much 
confidence  on  estimates  of  m * which  are  derived  from  small  populations  of 
data  points.  Because  of  the  small  number  of  constraints  available  in  such  cases, 
the  solution  will  not  be  greatly  overdetermined  and  will,  in  fact,  approach  a 
unique  solution. 

In  fact,  to  obtain  an  alternative  estimate  of  accuracy,  one  can  carry 
out  a simple  error  analysis  on  the  original  observation  equation  (2. 8).  In  this 
approach  the  errors  in  the  data  can  be  found  from  a knowledge  of  the  methods 
used  to  obtain  the  known  Ag's  (Brovar  et  al,  1964,  p.  278-279;  Kearsley,  1976, 
p.  120-122).  The  errors  in  the  covariances  could  be  estimated  by  statistical 
analysis  of  Ctj  in  the  course  of  the  evaluation  of  the  covariance  function  (Ibid, 
p.  98-100). 


2.4  Collocation  Theory 


In  collocation,  one  predicts  signals  from  observed  data  which  are  not 
necessarily  of  the  same  type  as  the  signal  predicted.  For  this,  auto-correlations 
and  cross-correlations  are  needed.  It  is  also  possible  to  incorporate  in  the  tech- 
nique the  determination  of  parameters  of  a geometrical  model  which  is  chosen 
a priori  in  order  to  reduce  systematic  effects  which  may  be  present  in  the  data. 

The  generalized  observation  equation  of  collocation  is  stated  as: 


(2.17)  x * AX  + s'  + n 

where 


x is  the  measurement  (e^.  gravity) 

AX  defines  a mathematical  model  (e.g.  normal  ellipsoid) 
s'  is  the  signal  (e.g.  gravity  anomaly),  and 
n is  die  noise  (e.  g.  measuring  error  in  gravity) 

Following  the  development  by  Moritz  (1972,  p.  7-16),  the  model 
parameters  X are  found  by  least  squares  to  be: 


(2.18) 


X - (AT  C’1  A)’1  At  C_l  x 
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and  the  predicted  signals  by: 


(2.19) 


s - Ctt  C"1  (x  - AX) 


where 


C - C,x,  the  auto-covariance  matrix  of  observations 

Csx  ■ the  cross-covariance  matrix  between  signal  and  observed  data . 


It  is  noted  that  the  expression  which  determines  the  model  parameters 
is  dependent  on  C and  thus  may  be  influenced  by  the  choice  of  stochastic  model . 

In  the  local  context  the  model  effects  are  generally  removed  by  adopting 
parameters  deduced  from  a "higher"  (or  global)  solution;  The  problem  then 
reduces  to  a "multivariate"  prediction  (Grafarend  and  Offermans,  1975,  p.  4), 
where  the  task  is  to  predict  one  potential-related  parameter  (e.g.  N or  £ ,77) 
from  a second  (e.g.  Ag).  It  is  this  situation  which  is  of  concern  in  this  present 
investigation. 

As  can  be  seen  from  (2. 19)  the  prediction  is  a function  of  both  the 
auto-covariances  of  the  observed  quantities,  and  the  cross -covariances  of  the 
observed  quantity  with  the  predicted  quantity.  Both  of  these  quantities  are  derived 
from  the  anomalous  potential,  and  their  covariance  functions  are  thus  indirectly 
related  (Moritz,  1972,  p.  94-99;  Grafarend  and  Offermans,  1975,  p.  28).  A 
number  of  models  have  been  suggested  as  suitable  for  representing  the  auto- 
and  cross-covariance  functions,  both  on  a global  scale  (Tscherning  and  Rapp, 

1974)  and  for  the  local  case  (Jordan,  1972).  However,  isotropy  of  the  potential 
field  and  hence  of  the  Ag  field  (Ibid,  p.  3664)  is  usually  assumed.  (See  also 
Grafarend  and  Offermans,  1975,  p.  20-22;  Grafarend,  1976,  p.  164-165). 

Thus  it  will  be  seen  that  the  values  resulting  from  (2. 19)  will  have  the 
characteristics  of  isotropy  imposed  on  diem  by  means  of  the  auto-  and  cross- 
covariance  functions.  It  is  possible  to  include  an  azimuth  term  in  the  functions 
' to  account  for  anisotropy  (Jordan,  1972,  p.  3663).  The  situation  is  much  more 
difficult  to  model  (see  section  4.3)  due  to  complications  which  occur  in  the 
cross-covariances  of  two  different  signals.  Nevertheless,  it  is  obvious  that  the 
adoption  of  an  isotropic  statistical  model  to  generate  covariances  in  a non-isotropic 
region,  while  simplifying  the  solution,  will  cause  smoothing  which  may  be 
detrimental  to  the  solution. 


2.  S Comments 


The  methods  of  least  squares  collocation  and  prediction  are  often 
described  as  the  "optimal"  solutions  of  their  respective  problems  given  the  data 
available  (Heiskanen  and  Moritz,  1967,  p.  269;  Moritz,  1972,  p.  18;  Lachapelle, 
1975a,  p.  19).  From  this  it  has  been  inferred  that  the  techniques  as  formulated 
provide  the  "best  possible"  method  of  predicting  or  solving  for  unknown  param- 
eters. Such  a claim  must  be  viewed  in  the  light  of  the  assumptions  made  in 
equations  (2. 2)  and  (2. 11).  If  a stochastic  model  with  properties  of  homogeneity 
and  isotropy  is  not  a good  description  of  the  actual  field  then  a solution  based 
upon  this  model  will  not  necessarily  give  the  most  accurate  solution,  hi  other 
words,  least  squares  prediction  or  collocation  will  give  the  optimal  solution  for 
the  particular  stochastic  model  adopted;  whether  or  not  this  is  the  best  of  all 
solutions  is  dependent  upon  the  faithfulness  of  this  model  to  the  real  situation. 

The  matter  as  to  whether  or  not  a more  sophisticated  covariance 
model  will  have  any  practical  impact  on  the  solution  has  already  been  touched 
on  in  section  1.  It  is  hoped  that  this  investigation  will  provide  some  further 
insights  into  this  question. 

3.  The  Detection  and  Representation  of  Anisotropic  Characteristics 
3. 1 Introduction 


In  this  chapter  various  methods  tested  as  a means  of  detecting 
anisotropy  are  reviewed.  Possible  ways  of  representing  the  azimuth-dependence 
of  the  statistical  relationships  are  also  suggested,  it  being  realized  that  the 
method  which  is  most  efficient  in  detecting  anisotropy  will  also  provide  the 
best  foundation  for  its  representation. 

To  aid  in  the  appreciation  of  the  methods,  tests  were  carried  out  on 
a number  of  data  sets,  details  of  which  appear  in  Table  1.  Because  of  the  limited 
extent  of  the  data,  it  was  sale  to  assume  plane  relationships  existed  between 
data  points  involved  in  the  analysis. 


Table  1 

[t  Details  of  Test  Data  Sets 


Data 
Set  No. 

Location 

Description 

Data 

1 

38°£^p£  39° 

-83°i  X*6  -84° 

Strong  ridge  (wavelength 
"“•30';  amplitude  »■ 

40  mgals) 
dominating  on  azi- 
muth of  20° 

Range  80mgal 
(see  Figure  1) 

Free-air  anomalies 
interpolated  from 
isogal  map  on  5' 
grid  interval 

2 

Generated 

Data 

Sloping  plane 

Max.  50mgal,SE  corner 
Min.  -2. 5mgal,  NW  cnr. 
Azimuth  of  strike*-  35° 

As  above 

3 

Generated 

Data 

Sloping  plane 

Max.  28mgal,  E boundary 
Min.  2mgal,  W boundary 
Azimuth  of  strike  0° 

As  above 

4 

38°*(p«  39° 
-82°a  Ai  -83° 

Flat,  featureless 

Range  40  mgal 

Note:  adjoins  data  set  1 
(See  Figure  2) 

As  above 

Special  attention  must  be  paid  to  the  dimensions  of  the  grid  in  any 
procedure  which  attempts  to  analyze  a continuous  field  represented  by  discrete 
points  at  grid  Intersections.  The  reader  is  referred  to  Blackman  and  Tukey 
(1959,  117-120),  Horton  et.al.  (1964,  590-591)  and  Nettleton  (1976,  159)  fora 
discussion  on  this  topic,  hi  this  particular  case  the  5’  grid  intersection  was 
chosen  as  giving  adequate  representation  in  an  already  smoothened  field.  In 
most  cases  the  maximum  'lag'  or  step  (r  in  equation  2.5)  was  limited  to  half 
the  number  of  data  points  in  each  row  or  column.  The  minimum  sample  size 
was  approximately  50.  lb  subsequent  analyses  the  r.  m.  s.  error  of  the  mean 
obtained  for  this  maximum  step  did  not  differ  significantly  from  that  for  the  initial 
step  sizes,  and  so  was  considered  a good  estimate  for  the  quantity  being  deter- 
mined here  (eg.  the  covariance,  coskewance,  etc.) 


| 1 

■ 

i;  m 

* t 


Figure  1.  Free-air  Gravity  Anomalies 
(Contour  Interval:  5 mgal) 


Data  Set  1 
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It  should  be  noted  that  these  data  sets  are  hilly  defined  - a situation 
which  will  rarely  occur  in  real  life.  Such  laboratory  conditions  are  needed  for 
this  investigation.  The  application  of  collocation  and  prediction  in  areas  with 
little  or  no  data  is  a separate  issue  and  is  beyond  the  scope  of  this  report.  It 
may  be  possible  to  use  cross-correlations  between  gravity  and  other  geophysical 
phenomena,  such  as  heat  flow,  to  extend  gravity  data.  For  further  information 
in  this  area  the  reader  is  referred  to  Kaula  (1967),  Woollard  and  Daugherty 
(1970;  1974),  Woollard  et  al.  (1975)  and  Groten  (1975). 


3.2  Covariance  Analysis  of  Profiles 


An  obvious  method  of  testing  a field  for  anisotropy  is  to  find  the 
covariance  functions  along  two  mutually  orthogonal  profiles  (or  ’transects', 
Whittle,  1954,  434)  of  the  field.  In  this  case  it  is  most  convenient  that  these 
profiles  be  taken  along  the  N-S  and  E-W  axes,  much  along  the  lines  adopted 
by  Rapp  (1964,  31-57).  Equation  (2.5)  is  therefore  modified  to  become: 


(3.1)  C(r)  - M [Agt  Agj  \ 

"CO  " CtT  i 

where 

U j * It  n 

“« ■ C} ,or  o 

and  Agi  is  centered  as  before  (see  Equation  2.7). 

The  functions  resulting  from  a profile  analysis  of  data  set  1 are 
illustrated  in  figure  3.  The  points  to  notice  are: 

(I)  the  steepness  of  the  E -W  function  cf.  the  N-S  fonctlon.  In 
fact,  C*  (r)  ■ 0 mgal3  at  r-  30  km,  while  at  this  separation 
Ch  (tf  *■  280  mgal3 ; also  G«  (r)  does  not  cross  C (r)  • 0 mgal3 
within  the  limits  of  the  analysis. 

(II)  the  location  of  the  general  covariance  analysis  (equation  2.5). 
This  does  not  take  up  the  mean  position  between  C*  and  Ct  as 
might  be  expected,  but  rather  tends  toward  (and  eventually 
crosses)  the  Ci  function.  This  phenomenon  will  receive  more 
attention  In  section  [(3.6)  (lv)]. 
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( ill)  Profile  analysis  taken  on  data  sets  2 and  3 produced  functions 
similar  to  those  for  data  set  1 (see  figure  4).  The  N-S  profile 
for  data  set  3 which  is  horizontal  produced  a horizontal  line 
reflecting  the  1:1  correlation  in  this  direction.  The  fact  that  . 
the  data  sets  1 and  2 also  display  greater  correlation  in  the 
N-S  than  in  the  E-W  direction  explains  why  their  CN  Sanctions 
decay  slower  than  do  the  Ct  curves. 

(iv)  One  sees  the  general  analysis  of  data  set  2 produces  a curve  which 
tends  to  follow  the  Ct  curve  as  for  the  first  data  set.  An 
explanation  for  this  will  be  given  below  (section  3). 


While  dramatically  illustrating  the  presence  of  anisotropy  this 
technique  does  not  help  to  quantify  this  property.  It  is  difficult  to  see  how  one 
may,  by  simple  adaption,  generate  elements  of  the  covariance  matrices  (unless, 
of  course,  the  data  points  lie  on  the  N-S  or  E-W  axis  passing  through  the  pre- 
diction point). 


3.3  Spectral  Analysis 


A spectral  analysis  of  the  data  can  provide  insight  into  the  wave-forms 
contributing  to  the  field.  A spectral  analysis  of  the  covariance  function  produces 
the  ’“power  spectrum"  of  the  field.  The  power  spectra  of  both  the  N-S  and  E-W 
profiles  may  provide  more  information  into  die  nature  of  the  anisotropy  of  the 
data  sets. 

For  the  general  background  on  this  technique  the  reader  is  referred 
to  Blackman  and  Tukey  (1957,  p.  53  and  121)  and  Miller  (1956,  p.  164-170). 
Applications  of  spectral  analysis  to  geodetic  or  geophysical  problems  can  be 
found  by  referring  to  Horton  et  al.  (1964),  Moritz  (1967),  and  Nettieton  (1976, 
p.  156-181).  A brief  outline  of  the  technique  will  be  given  here  to  assist 
understanding  of  the  resultant  spectra.  The  power  spectrum  S(oo ) is  the 
Fourier  Transform  of  the  covariance  Sanction: 


(3.2) 


C(T)  e 


t#T 


where  t is  the  lag  time  (the  equivalent  to  the  step  r in  equation  2. 5).  Now 
C ( r)  is  a symmetric  ftonction.  The  Imaginary  term  of  the  expanded  form  of  e~tU3T 
(i  sin  u;  T ) will  integrate  to  zero  on  multiplication  with  C (t)  to  produce: 
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C(  r ) [mGal*]  C ( r ) [m  Gal*]  C(  r ) [mGal*] 


FIGURE  4: COVARIANCE  FUNCTIONS  OF  PROFILES 


DATA  SET  2 


DATA  SET  3 


DATA  SET  4 
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=*  2 I C (T)  COS  i|)T  dT 

• 0 


This  can  be  farther  modified  as  C (t)  - 0 as  t-«, 
Hence,  for  all  practical  purposes: 

(3. 4)  S ( w ) = f C (T)  cos  u)  T d T 

* 0 


where  T is  the  period  of  the  covariance  function  C(T). 

When  the  covariance  function  is  represented  at  discrete,  regular 
intervals  (Ax),  equation  (3.4)  becomes  (Blackman  and  Tuckey,  1959,  p.  53 
and  121) : 


• - 1 

(3.5)  Sr  (cc)  = Ax  I^Co  + 2 Y C*  cos  p r rr/m  + C,cos  rrrj,  r = 0 - m 

p=  1 

where 

r is  the  harmonic  counter 

C0-C.  = C(0) - C(m),  and 

m is  the  maximum  lag  or  step  value 

and  Ax  is  the  lag  Interval,  in  this  case  1 grid  unit. 

It  is  sometimes  possible  to  get  meaningless  negative  values  for  the 
coefficients  of  S ( u:).  Smoothing  techniques  can  then  be  applied,  such  as  (i) 
"hanning"  (Ibid,  p.  53;  Davis,  1973,  p.  269): 

(3.6)  Sq1  ■ i (So  + Sj, ) 

S»l  * ^ St-tf  i St  + \ Sp+i  1 < p < m 

s.1  - i (S^i  + s,  > 


where  the  S1  coefficients  are  used  for  analysis;  or  (li)  replacing  the  coefficients 
Ct  in  (3. 5)  with: 

(3.7)  (0/2)  (1  + cos  (m/m)) 

(see  Horton  etal.,  1964,  p.  588). 
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The  power  spectra  for  the  N-S  and  E-W  profiles  of  data  sets  1,  2 and 
4 are  shown  in  Figure  5.  The  units  of  the  power  density  estimates  are  mgal 
squared  per  cycle  per  grid  unit.  The  abscissa  is  annotated  in  both  frequency (f- 
cycles  per  grid  unit)  and  harmonics  of  the  fundamental  wavelength  (shown  by 
arrows).  See  Appendix  Al. 

The  spectra  for  data  set  1 differ  dramatically  in  the  low  order 
harmonics,  the  spectrum  for  the  N-S  profile  (St)  maximizing  for  the  zero  order 
term  while  that  for  the  E-W  profile  (S<)  receives  its  maximum  contribution 
from  the  first  harmonic.  A similar  effect  is  noticed  in  the  spectra  for  data 
set  2.  While  this  clearly  shows  the  field  is  anisotropic  and  this  is  valuable 
information  to  have  for  (say)  the  removal  of  trend  surfaces,  it  is  not  useful 
in  the  evaluation  of  the  anisotropy  of  the  data.  In  fact,  these  spectra  must 
be  treated  cautiously  because,  due  to  the  limitations  placed  on  the  covariance 
analysis,  an  incomplete  wavelength  of  information  is  available  for  analysis. 

3.4  Semi- Vario grams 

The  se mi-va riogram  has  been  used  extensively  in  the  area  of 
interpretative  geology  (Matheron,  1963,  p.  1250;  1965).  Some  mention  has  been 
made  of  their  use  for  geodetic  purposes  (Monget,  1969;  Monget  and  Albuisson, 
1971),  but  their  application  in  this  field  has  been  limited.  This  function  is  closely 
related  to  the  structure  function,  used  by  VyskoEil  in  the  statistical  analysis  of 
gravity  in  Czechoslovakia  (VyskoSil,  1970). 

The  semi-variogram  is  defined  as: 

(3.8)  y(r)=  jJ^(Agl+r- AgOa  dS 


where  S is  the  area  under  consideration.  This  function  assumes  statistical 
properties  of  homogeneity  and  isotropy  and  its  similarity  to  the  covariance 
Amotion  is  noted.  It  is  in  fact  related  to  the  covariance  function  by  the  simple 
expression  (Ibid.,  p.  174;  Matheron,  1963,  p.  1253): 


(3.9)  y(r)  = C (0)  - C (r) 

This  relationship  helps  in  understanding  the  behavior  of  the  se  mi- 
va  riograms  for  profiles  of  the  data  sets  (Figure  6).  These  curves  are  the 
mirror  images  of  the  covariances  of  the  profiles,  originating  at  (0,0)  on  the 


S(co)  S(u>)  S(tu)-mgal2per  cycle  per  grid  unit 


Figure  5:  POWER  DENSITY  ESTIMATES  FROM  COVARIANCE  ANALYSIS 


DATA  SET  I 


DATA  SET  2 


* * i * 1 * | 
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DATA  SET  4 


0 2 3 4 5 6 
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r )/(mgal*) 


Figure  6:  SEMI-VARIOGRAMS  OF  PROFILES 


graph.  According  to  the  definitions,  they  fall  into  the  "continue -type"  class, 
l.e.  "represent  a regionalized  variable  of  high  continuity"  (Ibid. , p.  1250). 
While  again  demonstrating  the  presence  of  anisotropy,  they  do  not  provide  any 
information  not  already  known  from  the  covariance  functions.  In  fact,  the 
task  of  finding  elements  in  the  covariance  matrices  is  made  more  complicated. 
It  would  appear  that  their  main  value  lies  in  the  area  of  interpretative  geology. 


3.5  Coskewance  Functions 


The  coskewance  analysis  of  the  anomalous  gravity  field  was  suggested 
by  Kaula  (1966a;  1966b;  1967)  as  a method  of  performing  a non-linear  auto- 
analysis of  die  anomalous  field.  It  is  used  here  not  in  this  context  but  as  a 
means  of  detecting  and  measuring  the  anisotropic  nature  of  the  field.  Coskewance 
is  defined  as: 

(3.10)  L(T,u>,j9)  = L ((()„  , rt . 0r)  3 •— 7 JJJ x (r)  x (s)  x (t)  d t d s d r 


where  the  elements  are  shown  in  Figure  7. 


Figure  7 
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It  is  analogous  to  covariance  except  that  triple  products  rather  than,  product  pairs 
are  meaned  throughout  die  field. 

The  original  development  assumed  Isotropy  in  the  data,  i.e.  means 
were  taken  without  regard  for  Of  and  )3 . However,  if  ot  and  j9  are  kept  constant 
for  any  one  set  of  nbr«  and  ^rt  then  information  regarding  the  azimuth  dependence 
of  the  field  can  be  found.  For  example,  coskewances  can  be  found  for  N-S 
profiles,  where  a = 0,  /3  = 0;  0*  ^r,  s 0 s 4>rt  s and  E-W  profiles 
( O * 90;  0=  0;  0 * «ir,  * 444 ; 0 * tf)rt  s il> «•*)  ami  these  compared.  Alternatively, 

offset  configurations  can  be  taken  and  compared  eg.  Of  ® 0,  £ 3 90;  0 s 0r,  s ^«»x ; 

0s*rtS*,.,. 

In  this  study  the  parameter  4>r(  (or  its  planar  equivalent,  T)  is  called 
the  first  step,  while  tf>rt  is  the  second  step,  w.  An  example  of  the  curves  which 
result  from  the  coskewance  analysis  of  data  set  1 are  shown  in  Figure  8.  Note 
that  upon  assuming  isotropy,  L is  no  longer  a function  of  0. 


It  is  difficult  to  interpret  the  coskewance  curves,  although  the  effects 
of  anisotropy  are  apparent.  (If  the  field  were  isotropic  L (T»  u.')  would  equal 
L(ui,  t ) for  all  t,  u - a situation  which  clearly  does  not  exist  - particularly 
when  one  step  is  small  and  the  other  approaches  the  maximum.)  It  is  interesting 
to  note  that  L (3,  ui)  for  all  ware  approximately  equal,  reflecting  the  fact  that 
3 step  units  is  about  the  half  wave-length  of  the  data.  Also,  as  might  be  expected, 
L(T,  w)  for  w - 0 are  similar  in  shape  to  the  covariance  curve. 

However,  the  combination  of  three  parameters  in  this  fashion  does 
complicate  interpretation  and  gives  no  clear  indication  of  the  extent  of  anisotropy. 
It  is  also  difficult  to  use  these  functions  to  set  up  the  covariance  matrices  needed 
in  subsequent  computations. 


3. 6 Two  Dimensional  Covariance  Functions 


The  coskewance  functions  appeared  promising  because  it  combined  data 
in  a variety  of  configurations.  That  is,  it  involved  both  the  separation  and  azimuth 
as  parameters  in  the  analysis.  However,  the  addition  of  the  third  element  in  the 
products  distorted  the  desired  covariance  relationship  and  complicated  the 
interpretation  of  die  resultant  functions. 

With  a small  modification  - the  omission  of  pivot  point  r above  - one 
can  find  the  product  pairs  in  the  various  configurations  mentioned  in  section  3. 5. 
The  resultant  function  is  the  two-dimensional  (2-D)  covariance  function.  Thus: 

*-trl  w»|*l 

(3.11)  C(r,S)  - ■■  Y Y Ag(itj)  A«(i  + r,j+s),  r,s»0,m 

(N-|r|)(N-|s|)  1;1--1 
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where 


r,  s Is  the  step  in  the  row,  column  respectively 
m la  the  maximum  step  set  for  the  analysis  (usually  N/2),  and 
N is  the  number  of  rows,  columns  in  the  array  (assumed  square). 


This  expression  is  algebraically  correct,  that  is  the  steps  in  the  rows 
and  columns  can  be  both  positive  and  negative.  It  will  be  seen  that: 


C(-r,s)  * C (r,-s) 

C (-r,-s)  » C (r,s) 


In  the  resultant  surface  quadrants  I and  in  are  identical  as  are  quadrants  n and  IV. 

This  ftinction  has  been  found  useful  in  other  areas  of  scientific  study  to 
depict  the  azimuth-dependent  statistical  relationships  which  exist  between  data 
points.  For  example.  Whittle  (1954)  used  it  to  demonstrate  the  auto-correlation 
of  soil  fertility  as  evidenced  by  the  wheat  yield  of  rectangular  subdivisions  of  a 
test  area.  It  has  also  assisted  the  study  of  magnetic  anomalies  in  the  north- 
western regions  of  Canada  (Horton  et  al. , 1964,  p.  597-599).  It  has  been  used 
to  tBst  the  anisotropic  characteristics  of  the  gravity  field  in  the  Carpathian 
Mountains  of  Czechoslovakia  (Kubackova,  1974).  This  area  is  similar  in  size 
to  the  test  areas  in  this  present  study  and  data  was  taken  at  1'  x 1'  intersections. 
The  ftinction  has  received  mention  in  other  gravimetric  studies  (e.  g.  Jordan, 

1972,  p.  3661)  but  this  writer  has  found  no  other  use  of  it  in  the  geodetic  literature. 

The  2-D  covariance  analysis  of  data  set  1 produces  the  surface  Illustrated 
in  Figure  9.  The  following  points  should  be  noted: 


(i)  the  surface  represents  the  covariance  relationships  for  product 
pairs  of  all  separations  and  orientations.  Note  that  it  includes 
the  N-S  and  E-W  functions  deduced  above;  these  are  found 
by  taking  a section  through  the  midpoint  of  the  figure,  C(0,0), 
north  and  east  respectively. 

(11)  The  anisotropic  character  of  the  original  data  is  fully  quantified. 
The  predominant  direction  of  the  isogal  lines  in  figure  1 
(Qf=“  10°)  is  reflected  by  die  orientation  of  the  contours  of  the 
2-0  covariance  surface. 

(ill)  A truly  isotropic  data  set  will  produce  a circular  pattern  for  the 
contours  of  the  2-D  covariance  surface  (see  Agterberg,  1970,  125; 
Kuba&kova,  1974,  19).  The  extent  of  the  deviation  of  the  contours 
from  circles  indicates  the  extent  of  anisotropy  existing  in  the 
field. 
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(tv)  The  behavior  of  the  general  covariance  function  (section  3. 2)  can 
no«v  be  understood.  This  results,  of  course,  from  the  average 
of  all  profiles  0 * a<  360.  The  fact  that  the  contours  of  the  2-D 
covariance  surface  run  approximately  N-S  will  mean  that  the 
average  taken  for  any  one  separation  (r)  will  tend  toward  the 
curve  which  is  more  representative  of  the  predominant  trend  of 
the  surface  i.e.  the  Cc  curve. 

(v)  Because  of  population  differences  the  r.m.s.  of  the  C (L,  j) 
values  increase  with  increase  in  l,  j.  For  example,  the 
spread  for  C (0,0)  in  (the  r.m.s.  sense)  is  ± 30  mgal3  while 
the  same  spread  for  C (+6,  +6)  is  ±60  mgal3.  This  should  be 
remembered  when  extracting  C (i,  j)  values  for  subsequent 
computation. 

The  2-D  covariance  analysis  of  data  set  2 produces  a surface  with 
features  similar  to  the  surface  just  discussed  (see  Figure  10).  However, 
analysis  of  data  set  4 (adjacent  to  data  set  1)  shows  that  this  field  tends  much 
more  to  the  isotropic  situation  (see  Figure  11).  Insofar  as  it  bears  little 
resemblance  to  the  2-D  covariance  surface  for  data  set  1,  it  reinforces  the 
warning  concerning  the  use  of  the  covariance  function  derived  from  one  local 
field  to  represent  the  statistical  characteristics  of  another  field,  no  matter 
how  close  that  field  may  be. 


3.7  Representation  of  the  2-D  Covariance  Surface 


It  is  possible  to  represent  the  2-D  covariance  surface  in  three  ways: 
(i)  by  surface-fitting 
(11)  by  a family  of  curves 
(ill)  by  discrete  data  . 


Section  4. 2 shows  the  results  of  some  experiments  which  used  method 
(ill)  to  represent  die  2-D  covariance  surface  of  data  set  1 (Figure  9).  It  could 
be  well  represented  by  a fairly  low-order  surface  (say  of  order  3,  however, 
the  approach  which  uses  a family  of  curves  has  the  distinct  advantage  that  it 
directly  relates  to  the  functions  already  used  to  represent  general  covariances. 
This  advantage  is  enhanced  if  one  chooses  a function  which  is  compatible  for 
covariances  of  all  geoldal  relationships  (cf.  Jordan,  1972). 

The  basic  parameters  which  have  been  used  to  describe  the  normal 
(1-D)  covariance  Auction  are: 
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NORTHERLY  STEP 


NORTHERLY  STEP 


(1)  Co  - the  variance, 

(ii)  £ - the  correlation  length,  defined  such  that  C (0  = i Co,  and 
(111)  x - the  curvature  of  the  curve  at  r = 0 


(see  Moritz,  1976,  p.  22).  The  approach  developed  here  uses  Co  and  £ as  the 
basic  function  descriptors,  and  essentially  makes  the  £ parameter  a function  of 
azimuth.  In  this  way  it  is  possible  to  generate  a covariance  function  for  any 
azimuth  and  thus  derive  elements  of  the  covariance  matrix  needed  for  computations. 

The  techniques  used  in  two-dimensional  error  analysis  are  adapted  for 
this  purpose.  The  elliptical  shape  of  the  contours  of  the  (2-D)  covariance  surface 
has  already  been  noted;  it  is  possible  to  compare  the  curve  defined  by  £a  , 

0 < a s 360,  to  the  error  ellipse  of  error  analysis.  The  mathematical  techniques 
derived  for  the  transformation  of  the  error  ellipse  will  then  be  useful  in  our 
present  task. 

The  weight  coefficient  in  direction  Of  (Qaa)  can  be  expressed  as 
(Richardus,  1966,  p.  100-105;  Hirvonen,  1971,  p.  165-169): 


(3.12)  Qaa  = Qxx  cos3  a + Q,y  sin3  at  + 2QXy  sin  a cos  a 


Qxx » Qyy  = weight  coefficient  along  the  X,  Y axis  respectively 
Q*y  ■ the  cofactor  between  X Y coefficients 


Now,  defining  £N  , £t  as  the  £ value  for  C* , Ce , and  substituting  them 
for  Qxx  , Qyy  in  (3.12): 


(3.13)  £q  = £n  cos3  o + £i  sin3  Of  + 2 £N£  sin  a cos  a 


The  value  £ni  needs  to  be  defined.  Being  a cofactor  there  is  no  geo- 
metrical meaning  attached  to  it  as  there  is  for  £N  , £c  • However,  it  can  be 
determined  if  the  correlation  length  is  known  for  at  least  one  azimuth  other  than 
0°  or  90°.  Then  by  solution  of: 

(3.14)  £nc  » UoraM^N  cos3  a + £f  sin3  a)]/2  sin  a cos  a 
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this  quantity  can  be  determined.  The  most  sensitive  estimate  would  be  obtained 
for  a mid-quadrant  azimuth,  and  to  control  the  discrepancies  occurring  due  to  the 
departure  of  the  £ curve  from  a true  ellipse,  two  values  of  £Nt  were  obtained 
from  a * 45°  and  a * 135°.  The  mean  of  these  was  then  substituted  into  (3. 13). 
This  expression  could  then  be  used  to  generate  the  covariance  function  for  any 
azimuth. 


Local  Covariance  Functions 


Many  expressions  have  been  suggested  to  represent  the  local  covariance 
flinction.  Experiments  were  made  on  a selection  of  these  to  determine  which  was 
best  suited  for  surface  generation.  Those  tested  are  listed  below.  The  parameter 
s represents  the  step  or  lag  throughout. 

Model  1:  Hirvonen’s  covariance  function  (Hirvonen,  1962;  Moritz,  1976,  p.  17) 


(3.15) 


C (s)  = C«/{l  +(s/d)^ 


For  this  function  £ = d 

Model  2;  The  Gaussian  function  (Ibid. , p.  18,  26) 
(3.16)  C (s)  = Co 


For  which  £ = 


Models  3 and  4; 


These  were  Models  1 and  2 compounded  with  the  cos  term,  cos  (rrr/A). 
This  was  done  in  order  to  obtain  the  characteristic  negative  values  for  the  covar- 
iance function.  Unmodified,  Models  1 and  2 approach  zero,  but  never  become 
negative.  The  wavelength  (A)  of  6 units  as  used  in  the  spectral  analysis  was 
again  adopted.  However,  difficulty  was  encountered  in  maintaining  the  £ values 
for  this  wavelength,  and  the  results  obtained  were  generally  poor.  They  are 
omitted  from  the  summarized  results  in  Table  2. 
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Model  5;  The  logarithmic  model  (Moritz,  1976,  p.  29) 


(3.17) 


C (a)  * fti 


2 e 


A l+(l+kass)^ 


where  £ - k-1  <[(2eA/^-l)3  -l}1 


Assuming  k * 1 to  maintain  the  accepted  values  for  the  curvature 
parameter  (Ibid.). 


A = 2 1 


*1 


Model  6;  The  third-order  Maikov  model  (Jordan,  1972,  p.  3664) 
(3.18)  C(s)  = Co  (1  + ^ -55a)  e-/o 


£ is  found  to  be  approximately  1.095D  (Ibid. , p.  3665).  This  function 
is  of  interest  because  of  its  compatibility  with  the  auto-  and  cross-correlation 
functions  of  other  derivatives  of  the  anomalous  potential. 

Using  all  models,  C(s)  was  found  for  all  s,  l*ss5,  with  £a  computed  from  (3. 13) 
for  azimuths  0 s a < 360  in  10°  advances.  A selection  of  these  are  tabulated  in 
Table  2,  and  compared  with  the  actual  discrete  values  of  the  2-D  covariance. 

Also,  out  of  interest,  the  values  for  the  general  covariance  fiinction  are  listed. 

The  comparisons  indicate  that  the  logarithmic  model  (model  5)  appears 
to  give  the  best  all  around  results.  All  models  fit  well  for  a = 0 where  the  curve 
attenuates  slowly.  The  poorest  fit  is  found  for  negative  covariances,  and  it  is  in 
this  region  that  model  5 is  clearly  superior.  It  must  be  noted  that  discrepancies 
In  azimuths  50°  and  130°  will  be  influenced  by  the  approximation  of  the  $ - generator 
(equation  3. 13)  to  the  actual  £ curve.  Tbeat  discrepancies  could  be  minimized 
by  the  adoption  of  a more  sophisticated  expression  to  determine  £a. 

It  appears  that  this  technique  is  quite  successful.  It  must  be  pointed 
out  that  generation  of  covariances  in  this  way  comes  much  closer  to  the  truth 
than  the  values  found  from  the  general  covariance  fiinction  based  on  an  isotropic 
statistical  model.  Whether  or  not  this  improvement  in  accuracy  is  significant 
must  be  the  subject  of  a further  investigation  (section  4. 1). 


32 


3.8  Conclusions 


The  2-D  covariance  Amotion  promises  to  be  the  most  efficient  way  of 
detecting  and  representing  anisotropy  in  a local  or  regional  data  set.  It  clearly 
illustrates  the  extent  to  which  the  data  departs  from  an  Isotropic  situation. 
Because  it  represents  the  covariance  relationships  which  exist  between  point 
pairs  of  all  separations  and  orientations,  it  is  easily  used  to  generate  elements 
in  the  covariance  matrices  for  prediction  or  collocation  computations. 

The  2-D  covariance  surface  can  be  represented  by  means  of  a family 
of  curves,  generated  using  Co  computed  in  the  analysis,  and  the  correlation 
length  ^ found  for  four  azimuths  from  the  covariance  surface.  The  logarithmic 
model  of  Moritz  (model  5)  appears  to  represent  the  surface  most  successfully. 

It  is  thought  that  an  improvement  in  the  £ generator  would  similarly  improve 
the  covariance  values  in  the  mid-quadrant  directions. 


4.  The  2-D  Covariance  Function  in  Collocation  and  Prediction 


4. 1 Introduction 


The  2-D  covariance  function  represents  the  covariance  relationships 
of  point  pairs  of  all  separations  and  orientations.  The  question  now  arises  as 
to  the  possibility  of  using  this  function  in  collocation  and  prediction  computations. 
Most  theoretical  developments  which  derive  analytical  expressions  to  account 
for  an  anisotropic  situation  modify  these  by  assuming  Isotropy  before  continuing 
with  the  investigation  (e.  g.  Jordan,  1972,  p.  3663;  Grafkrend  and  Offermans, 

1975,  p.  14).  This  occurs  even  for  the  cross-covariances  of  the  deflections  of 
the  vertical.  These  are  shown  to  be  azimuth  dependent  whether  or  not  the 
potential  field  is  assumed  isotropic  (Moritz,  1972,  p.  Ill;  Grafarend  and 
Offermans,  1975,  p.  14).  However,  it  should  be  noted  that  if  the  potential  field 
is  assumed  isotropic,  the  expressions  for  the  above  cross -covariances  are,  in 
fact,  a modification  of  the  more  complete  expressions  which  account  for 
anisotropy  (Ibid. , p.  15). 

However,  there  appears  to  be  no  theoretical  objection  to  the  use  of 
the  2-D  covariance  function  in  collocation  or  prediction  solutions.  The  compu- 
tation is  thought  to  be  in  a planar  rather  than  a linear  (or  time)  domain.  The 
contribution  of  the  known  element  to  the  unknown  is  now  dependent  on  its 
directional,  as  well  as  its 'temporal',  position  on  the  plane.  The  apt  in 
equation  (2.14)  accounts  for  this,  for  the  elements  of  the  array  Cpi  are  now 
derived  having  regard  to  the  relative  position  of  the  signal  point  vis-a-vis  the 
known  point.  Similarly,  the  CtJ  matrix  is  developed  by  noting  the  disposition 
of  each  data  point  to  all  other  data  points  in  turn.  These  values  can  be  generated 
using  techniques  described  in  section  3. 8.  However,  in  the  experiments  below, 
values  have  been  taken  directly  from  the  arrays  resulting  from  the  2-D  covariance 
analysis. 
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4. 2 Re aiil ta  of  Prediction  Testa 


The  aim  of  the  teste  was  to  compare  results  of  predictions  made  using 
the  2-D  covariance  function  with  those  resulting  from  use  of  the  general  covariance 
function.  Data  sets  1 and  2 were  used  in  these  tests  as  these  were  two  fields 
exhibiting  strong  azimuth  dependent  characteristics. 

In  tbs  test  eight  "known”  points  were  used  to  estimate  an  unknown  point. 
(This  number  if  close  to  the  optimum  suggested  by  Rapp  (1964,  p.  141)  for  the 
prediction  of  point  anomalies.)  These  eight  points  were  chosen  in  various 
configurations  in  order  to  exaggerate  the  effect  of  the  anisotropy  of  the  field  (see 
Figure  12).  This  configuration  was  "stepped"  through  the  data  array,  and  the 
differences  between  true  and  predicted  values  used  to  find  the  actual  r.  m.  s.  error  of 
prediction.  This  value  was  also  compared  with  the  theoretical  r.m.s.  error  derived 
from  equation  (2. 16).  See  Table  3 for  definitions. 


4.2.1  Configuration! 
(1)  Data  Set  l 


The  results  of  prediction  using  the  2-0  covariance  function  are 
marginally  worse  than  those  using  the  general  function  (see  Table  3).  The 
distribution  of  errors  is  similar  (see  Figure  12), although  an  inspection  of  a 
plot  of  these  differences  shows  that  the  2-D  function  consistently  produces 
larger  discrepancies.  It  appears  that  the  data  points  closest  to  P (the  pre- 
diction point)  have  the  greatest  influence  on  the  predicted  value.  Their  prox- 
imity to  p (1. 4 step  units)  meant  that  there  was  no  significant  difference  between 
the ’weights'  produced  from  the  2-D  and  the  general  covariance  functions.  For 
example,  Cm  (1.4)  * 260  mgal*  of  Cm  (1,1)  3 245  nagal3  and  Cm  (1,-1)  3 240 
mgal.  It  could  be  argued  that  a lack  of  homogeneity  throughout  the  field  meant 
that  the  2-D  covariance  function  was  unrepresentative  of  the  average  situation 
throughout  the  whole  data  set,  and  it  was  adequate  to  assume  a hilly  random 
statistical  relationship  between  point  pairs  at  this  small  separation. 

(11)  Data  Set  2 

In  d»i«  case  the  effect  of  anisotropy  is  felt.  The  2-D  function 
produces  an  r.  m.  s.  of  ± 11. 6 mgal  (and  m,  “ ± 7. 1 mgal)  while  die  general 
hmetion  gives  rise  to  an  r.m.s.  of  ±13.4  mgal  (and  an  mp  of  ±11.7  mgal). 

This  is  not  unexpected  in  this  extreme  case.  The  fact  that  the  difference  is 
not  larger  must  in  part  be  due  to  the  proximity  of  the  closest  points  to  P. 

(See  also  the  distribution  of  differences.  Figure  12). 
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FIGURE  12  : RESULTS  OF  PREDICTION  TESTS  [umts:mgai] 
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4.2.2  Configuration  2 


LI)  Data  Set  1 


Configuration  2 was  deliberately  chosen  to  maximize  the  difference 
between  the  general  function  and  the  2-D  function.  For  example,  C (3, 1)  = 240 
mgal3  and  C (4,0)  - 220  mgal3  while  the  general  function  equivalents,  C (3.2), 
C (4),  were  70  mgal3  and  -12  mgal3  respectively.  Results  of  the  tests  showed 
that  the  anisotropic  model  fitted  the  data  better  (mp  * * 7 mgal  for  the  2-D 
model  and  ± 12  mgal  for  the  isotropic  model).  A similar  Improvement  existed 
in  the  actual  error  of  prediction  (±14  mgal  cf.  ±16. 5 mgal).  This  supports 
the  contention  that  the  choice  of  an  isotropic  model  does  not  always  produce  the 
optimal  solution.  It  also  underlines  die  danger  of  using  the  parameter  mp 
without  reservation  as  an  estimate  of  the  error  of  prediction. 

(ii)  Data  Set  2 


The  2-D  covariance  function  produces  a very  low  'actual'  error 
(±  0.6  mgal),  mp  being  computed  at  ±4  mgal.  However  the  solution  using  a 
general  covariance  function  breaks  down.  Individual  errors  of  up  to  100  mgal 
are  produced  and  the  value  for  mp  shows  a large  error  of  prediction. 
Apparently,  it  is  not  possible  to  perform  any  meaningful  prediction  (with  an 
asymmetric  configuration)  using  the  general  function,  when  conditions  such  as 
those  in  data  set  2 exist. 


Table  3 


Statistical  Analysis  of  Prediction  Tests 


Configuration  1 j Configuration  2 


Data  Set  1 Data  Set  2 I Data  Set  1 Data  Set  2 


General  2-D  General  I 2-D  General  2-D  General 


11.6  13.4  I 14.0  16.5  0.6  7.1(?) 


11.3(7) 


6.2  4.6  5.8 


1.16  1.16  1.05  0.68  I 1.07  1.09  1.04  -5.48 


(*)  - (predicted  Ag  - actual  Ag)a/no.  of  predictions  ? 
(+)  - error  of  prediction  - see  equation  (2.16) 
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Attention  ia  now  directed  toward  the  third  row  of  Table  3 where  the 

sums  of  die  coefficients: 
t* 

i 


a*i  * C*I  Ci? 


(derived  from  equation  2. 14)  are  listed.  These  coefficients  are  basic  in  the 
evaluation  of  both  the  predicted  value  of  Ag  (equation  2. 8)  and  the  theoretical 
estimate  of  accuracy  (equation  2. 16).  It  is  of  interest  to  note  therefore,  that 
in  general,  the  most  accurate  results  were  obtained  from  that  covariance  function 

giving  the  ) a»t  closest  to  unity. 

** 

1 


It  is  difficult  to  know  the  actual  significance  of  this.  Equation  2. 8 
looks  similar  to  the  normal  "weighted  mean"  in  this  context,  but  is  not  identical 
to  it.  It  has  already  been  noticed  that  V a,t  is  not  necessarily  unity  (Rapp,  1964, 
p.  7).  In  some  cases  it  is  valid  to  adjiist  the  individual  elements  of  the  an  vector 
to  let  the  sum  equal  unity.  In  this  regard  it  is  worth  quoting  from  Holloway  (1958, 
p.  354-355)  "The  sum  of  the  weights  of  a smoothing  or  filtering  fimction  determines 
the  ratio  of  the  mean  of  the  smoothed  or  filtered  series  to  the  mean  of  the  original 
series.  In  smoothing  it  is  generally  desired  to  leave  the  mean  of  the  series 
unchanged,  so  ^ a?{  made  equal  to  unify.  With  some  filters  it  is  not  necessary 

to  preserve  the  mean  of  the  series  and  in  these  cases,  ^ a»t  # 1. " 

i 

hi  the  prediction  context  it  is  desirable  that  the  property  of  ergodtelty 
hold.  It  therefore  appears  that  the  covariance  function  which  produces  the  Y a i 

i 

nearest  unity  is  the  one  which  maintains  this  property  best;  hence  it  produces  the 
best  comparisons  of  predicted  values  versus  known.  On  the  other  hand,  computations 
with  ^a,  made  equal  to  unity  showed  that  some  caution  must  be  observed  when 

i 

modifying  the  coefficients  in  this  fashion.  This  was  particularly  the  case  when 
the  Y a’i  differed  greatly  from  unity  (e.g.  0.50).  In  this  case,  although  the 
» 

bp  value  Improved,  some  predictions  changed  markedly  ("“  30  mgal)  and  the 
calculated  r.m.s.  deteriorated.  (See  starred  items,  Table  4,  Section  4.3.) 

However,  when  the  difference  between  £ a»t  and  unity  was  small  and  such 

i 

modifications  made,  a slight  Improvement  in  prediction  was  noticed  (again  see 
Table  4). 

I 
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4. 3 Results  of  Collocation  Tests 
4.3.1  Introductory  Remarks 


The  area  chosen  for  the  tests  was  the  region  near  the  SC  coast  of  the 
United  States  used  for  calibration  of  Geos-3.  It  is  bounded  by  40  * <p  > 20  and 
297  2X2  277,  the  information  on  latitude  20°  being  deleted  because  the  larger 
anomaly  values  here  introduced  erratic  behaviour  into  the  covariance  functions. 

The  gravity  data  was  in  the  form  of  1°  x 1°  mean  free-air  anomalies  extracted 
from  die  file  of  world-wide  1°  means  held  at  the  Department  of  Geodetic  Sc  ience. 
The  separations  (N)  have  been  computed  in  this  area  using  the  GEM  6 potential 
coefficients  to  degree  16  and  l°x  1°  gravity  data  (Rapp  and  Rummel,  1975).  Repre- 
sentations of  these  data  sets  are  illustrated  in  Figures  13  and  13a,  respectively. 

The  aim  of  the  collocation  test  was  to  compute  the  geoidal  separation 
(N)  from  the  gravity  anomalies.  The  computation  can  be  thought  of  as  a "multi- 
variate" prediction  as  defined  by  Grafarend  (1976,  p.  152)  because  the  model 
effects  are  already  assumed  to  be  removed.  (In  other  words  the  parameters 
AX  in  equation  2. 19  have  been  computed  independently  and  die  effects  of  normal 
gravity  removed  from  the  observed  gravity  to  obtain  the  free-air  anomalies.) 

The  computations  must  be  considered  preliminary,  although  the  author 
considers  them  to  be  a good  approximation  to  a more  rigorous  solution.  The 
reasons  for  this  are  as  follows: 

(1)  Sphericity  of  the  earth  (i.e.  die  convergence  of  the  meridians) 
was  not  accounted  for  and  the  data  was  assumed  planAr. 
Computations  of  the  general  covariance  and  profile  covariance 
functions  assuming  a spherical  model  showed  no  apprec  able 
change  from  the  same  functions  derived  using  a plane  model. 

Any  difference  was  certainly  minor  when  compared  with  the 
difference  found  between  the  general  covariance  functions  and 
the  2-D  functions. 

(li)  The  mean  anomalies  were  treated  as  if  they  were  point  anomalies 
In  the  analysis  process,  as  described  by  Smith  (1974,  p.  33).  All 
auto  and  cross-covariances  were  computed  from  the  known  data. 
That  is,  no  attempt  was  made  to  represent  the  covariances  from 
previously  determined  models. 

(iii)  It  should  be  remembered  that  this  test  is  to  obtain  a comparison 
of  the  2-D  function  against  the  1-D  function.  Any  shortcomings 
due  to  the  approximations  mentioned  above  will  affect  both  equal- 
ly. The  validity  of  the  comparison  will  therefore  not  be  harmed 
to  any  significant  degree. 
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4. 3. 2 Auto  and  Cross -Covariance  Analysis 


The  «u to -co variance  functions  resulting  from  a 2-D  analysis  of  the 
gravity  (data  set  5)  and  undulations  (data  set  6)  in  the  test  area  are  shown  in 
Figures  14  and  15.  (It  Is  interesting  to  note  that  the  strong  anisotropic  effect 
in  the  gravity  field  is  reflected  in  both  of  these  functions. ) The  contours 
below  200  mgala  show  no  signs  of  converging  within  the  limits  of  the  analysis. 
(At  step  10,  Cn  Is  +140  mgal3 . Ci  * 0 at  450  km.)  This  results  from  the 
obvious  N-S  trend  in  the  data  and  is  comparable  to  the  effect  seen  in  the 
analysis  of  die  N-S  profile  of  data  set  m (Figure  4).  The  angular  nature  of  the 
contours  must  be  due  to  the  use  of  the  1°  means  without  applying  any  smoothing 
to  this  data. 


The  2-D  covariance  surface  for  N (Figure  15)  indicates  the  presence 
of  anisotropy  in  this  data  and  again  reflects  the  orientation  of  the  trend  in  the 
original  data  (Figure  13). 


The  results  of  the  cross -covariance  analysis  are  shown  in  Figure  16. 
This  is  found  by  amending  equation  3. 11  to  include  the  two  variates.  Thus: 

N—  r h-  ■ 


(4.1) 


Cm.U  (*.s) 


1 

(N-|r|)(N-|s|) 


X £ N(l»J)  AS(i+r.  j+s) 

t»i  j=i 


when  applied  to  data  sets  5 and  6 produces  the  surface  represented  by  Figure  16. 
The  strong  anisotropic  trend  is  still  obvious,  although  the  angularities  of  Figure  14 
have  largely  disappeared  in  this  combination  of  gravity  data  with  undulations. 

Two  other  pecu’arities  are  noticed. 

(i)  The  ftinction  is  no  longer  diametrically  symmetric,  i.e.  quadrant  1 
i quadrant  3 and  quadrant  2 i quadrant  4.  The  reason  for  this  is 
apparent  on  looking  at  equation  (4. 1).  The  mean  of  cross-products 
of  N with  Ag  in  any  one  direction  does  not  equal  (except  accidentally) 
the  mean  of  cross-products  of  Ag  with  N in  that  direction.  If  the 
mean  of  these  two  quantities  is  used,  die  surface  resulting  will  be 
symmetric  (Figure  17).  However,  to  be  strictly  correct  the  first 
analysis  should  be  used  in  collocation  computations. 

(ii)  The  function  does  not  necessarily  behave  in  the  same  way  as  the 
auto-covariance  Auction  insofar  as  C(0,0)  is  not  necessarily  the 
maximum  value  achieved  by  the  surface.  For  example,  whereas 
C(0,0)  « 101.4  *9.8,  C(1,0)  - 106.1  * 9.5;  C(0,1)  « 103.4*9.2 
and  C(-2, 1)  * 104. 2*9.0.  The  author  has  seen  no  mention  of  this 
phenomenon  in  the  literature  and  it  is  something  which  should  be 
recognized  as  a possible  characteristic  of  the  cross-covariance 
analysis.  However,  because  the  increase  from  C(0,0)  is  small 
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Figure  14:  2-D  Auto-covariance  Function  - Data  Set  5 
(Contour  interval:  100  mgala) 
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Figaro  15:  2-D  Auto-covariance  Function  for  N 
(Contour  interval:  20ma) 
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Figure  17:  Mean  C roes -covariance  Function  Between  N and  Ag 
(Contour  interval:  50  mgal  m) 
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(particularly  when  viewed  against  the  standard  errors  of  the 
covariances)  it  seems  reasonable  to  assume  a small  decay  in 
the  ftxnction  to  assist  modelling  of  the  surface. 


Jordan  (1972,  p.  3665)  has  developed  expressions  which  relate  the 
variances  and  correlation  lengths  of  the  covariance  functions  for  the  gravity 
and  undulation  data.  It  would  be  ideal  if  such  relationships  could  be  used  to 
generate  the  2-D  covariance  functions  of  all  geoidal  parameters  from  (say)  the 
function  for  Ag  using  the  technique  described  in  section  3. 7.  Unfortunately, 
preliminary  calculations  showed  that  these  theoretical  relationships  did  not 
describe  the  actual  situation  at  all  well.  It  may  be  that  this  data  is  not  localized 
enough  for  such  modelling  and  that  it  is  necessary  to  use  functions  more  global 
in  nature  (e.g.  T scheming  and  Rapp,  1974). 

The  modelling  of  2-D  auto  and  cross-covariance  functions  needs 
further  investigation.  However,  at  this  stage  it  was  considered  more  important 
to  discover  whether  the  use  of  the  2-D  functions  had  any  significant  effect  in 
collocation.  For  this  reason,  discrete  data  resulting  directly  from  the  analysis 
of  the  local  field  was  used. 


4. 3. 3 Results  of  Tests 

a)  Univariate  Predictions:  Ag  from  Ag 

For  interest,  a simple  prediction  of  Ag  using  the  2-D  auto-covariance 
function  was  performed,  and  the  results  compared  with  a similar  computation 
using  the  general  covariance  analysis  (not  illustrated).  The  results  are 
summarized  in  Table  4 below.  Configuration  2 refers  to  the  point  arrangement 
described  in  Figure  12,  and  Configuration  3 is  described  in  Figure  17. 


Table  4 


Summary  of  Results  of  Predictions  - Data  Set  S 
(units  - mgal3) 


Covariance 

Function 

ESHEfflBB 

population  = 156) 

2-D 

General 

2-D 

General 

2-D 

General 

warn 

General 

lTe  XXI*  3* 

error 

15.1 

15.2 

15.1 

15.1 

16.3 

15.6 

21.7* 

15.4 

m 

E£n 

22.3 

19.3 

21.4  1119.2 

22.1 

13.0 

21.5 

mm 

HI 

mm 

+1.0 

Ha 

H9 

0.79 

HH 

H 

t Y a » forced  to  1. 0 


t 

* see  text  for  comment 
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Note  that  the  last  two  columns  (or  both  configurations  list  the  results  of  computations 

after  Y a t has  been  adjusted  to  unity.  Some  improvement  is  noticeable  generally 

in  both  the  r.  m.  s.  error  and  m> . The  exception  is  the  r.  m.  s.  error  for  the  2-D 
function  in  configuration  3 (marked  *).  Here  the  r.  m.  s.  deteriorates  and  individual 
differences  are  found  to  Increase  by  as  much  as  30  mgal. 

The  results  are  in  accord  with  those  above  (section  4.2).  Configuration 
3 is  symmetrical  about  the  prediction  point  and  the  general  function  produces  the 
superior  accuracy.  Larger  differences  occur  in  the  2-D  case  where  the  data  is 
not  well  represented  by  the  2-D  ftxnction  (particularly  in  the  SE  corner)  and  this 
is  a big  factor  in  producing  the  inferior  r.  m.  s.  error  overall.  If  the  prediction 
was  limited  to  the  regions  which  displayed  the  general  trend,  it  appears  that 
because  the  differences  are  smaller,  the  2-D  covariance  function  would  produce 
the  superior  results. 

The  2-D  and  general  (Unction  appear  to  be  equally  accurate  in  the  case 
of  Configuration  2.  The  population  for  this  configuration  is  large  (<»230)  and  15 
mgal  seems  to  be  the  limit  of  accuracy  for  this  computation.  It  is  noticed  that 
the  2-D  (Unction  converges  more  quickly  to  this  limit,  suggesting  (as  does  its  m* 
value)  a superior  accuracy  for  this  covariance  (Unction. 


Multivariate  Prediction:  N from 


The  value  of  N was  computed  throughout  the  field  from  the  gravity 
data  by  use  of  equation  2. 19,  with  parameters  AX  set  equal  to  zero.  The  gravity 
points  were  organized  in  Configuration  2 for  the  first  test.  Configuration  3 for 
the  second.  The  cross-covariances  used  to  provide  the  elements  of  matrix  Csx 
were: 

(i)  the  strict  2-D  cross -covariance  function  shown  in  Figure  16 

(ii)  the  averaged  or  symmetric  2-D  (Unction  (Figure  17),  and 

(Lli)  the  general  cross-covariance  (Unction  (not  illustrated) 

This  produced  three  solutions  for  each  of  die  configurations.  The 
elements  of  the  C matrix  were  taken  from  Figure  14  for  solutions  (1)  and  (li), 
and  from  the  general  auto-covariance  analysis  of  die  gravity  for  solution  (ill). 

A histogram  of  the  absolute  differences  |v|  resulting  for  each 
solution  is  shown  in  Figure  18,  and  the  comparisons  of  the  accuracies  achieved 
are  tabulated  in  Table  5.  Note  that  v is  found  by  differencing  the  predicted 
N value  from  the  N value  obtained  from  Figure  13a. 
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Table  5 


1 Summary  of  Accuracy  of  Collocation  Computations 

I (Units  : Meters) 


Configuration  2 

■a 

2-D 

Symmetric 

General 

■a 

2-D 

Symmetric 

General 

r.m.  s.  ± 
error 

4.71 

4.74 

4.83 

4.06 

4.00 

3.96 

rnr  ± 

8.0 

8.0 

8.1 

7.6 

7.8 

7.7 

I*. 

1 

0.39 

0.39 

0.45 

0.49 

0.48 

0.53 

The  comparison  of  predicted  versus  known  values  using  Configuration  2 
shows  that  the  strict  2-D  covariance  analysis  gives  the  better  results,  producing 
an  r.m.s.  of  4 4.71  m as  compared  with  ± 4.83  m for  the  general  case.  (It  should 
be  remembered  that  the  units  of  the  predicted  quantity  are  meters  and  a difference 
of  0.1  in  the  r.m.s.  is  not  insignificant.)  This  situation  is  reflected  also  by  the 
nw  values  (±8.0  met  4 8.1  m for  die  general  case). 

For  Configuration  3 the  general  covariances  produce  the  superior  results 
(±3.96  m as  against  ± 4.06  m for  the  2-D  covariance  values).  However,  the 
theoretical  accuracies  do  not  agree  with  this.  This  suggests  that  although  the 
2-D  statistical  model  is  a better  descriptor  of  the  field  overall,  larger  errors 
are  produced  by  this  model  in  uncharacteristic  areas  of  the  field  (see  Figure  18). 
Analysis  of  the  distribution  of  differences  from  both  the  general  and  2-D  solutions 
shows  a similar  pattern  exists  for  both,  and  that  large  errors  (|v|  > 5 m)  tend  to 
occur  where  the  gravity  profile  through  the  computation  point  increases  markedly 
away  from  the  point.  Neither  the  2-D  nor  the  general  covariance  analyses  are 
good  reflections  of  this  situation.  However,  the  general  function  seems  to 
exaggerate  this  effect  less  and  produces  the  smaller  difference. 

However,  in  areas  where  the  trend  of  the  field  is  strong  (e.  g.  around 
the  middle  western  area  of  the  field)  the  2-D  ftinction  produces  consistently  better 
results  from  both  configurations.  Here  anisotropy  is  present  and  it  would  be 
expected  that  the  2-D  analysis  prove  superior. 

It  is  also  of  interest  to  note  that  the  symmetrical  2-D  covariance 
function  produces  comparisons  between  those  obtained  by  the  strict  2-D  and  the 
general  ftinction,  as  may  be  expected.  The  summation  of  the  coefficients  (row  3, 
Table  S)  are  related  to  the  values  of  m»  (as  they  must  be),  but  the  significance 
of  the  value  (“0. 4)  is  not  clear.  Generalizing  from  earlier  comments  in  Section 
4. 2,  this  value  must  relate  to  the  ratios  of  the  mean  values  of  the  two  data  domains 
i.  e.  the  known  and  the  predicted.  The  ratio  of  the  mean  value  of  gravity  to  the 
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mean  value  of  N In  the  test  area  is,  in  fact,  0.46.  However,  die  relationship  noted 
in  the  prediction  computations  between  the  £ &t  value  and  the  accuracy  does  not 

appear  to  hold  here.  That  is,  the  greater  accuracy  does  not  necessarily  come  from 
that  computation  whose  sum  is  closer  to  0. 46. 

There  is  obviously  not  enough  information  here  to  come  to  any  definite 
conclusions  concerning  the  use  of  the  2-D  function  in  collocation.  Nevertheless, 
it  does  seem  to  follow  the  pattern  noticed  from  the  prediction  results.  That  is, 
where  the  anisotropy  portrayed  by  the  2-D  covariance  analysis  is  present,  the 
values  computed  from  the  2-D  covariances  compare  more  favorably  with  the  known 
data  than  do  those  computed  from  the  general  analysis. 


5.0  Conclusions 


The  two-dimensional  covariance  ftmction  provides  the  most  efficient 
means  of  detecting  and  representing  the  anisotropic  characteristics  of  a dat*  set 
distributed  over  a plane.  This  ftmction  graphically  describes  the  covariance 
which  exists  between  point  pairs  of  all  separations  and  orientations.  The  extent 
of  anisotropy  is  indicated  by  the  departure  of  the  contours  of  the  2-0  covariance 
surface  from  a circular  pattern,  and  the  orientation  of  the  axes  of  maximum  and 
minimum  correlation  are  clearly  shown. 

It  is  possible  to  model  the  2-0  covariance  surface  by  generating  a simple 
covariance  ftmction  for  each  azimuth,  0 * a < 360.  The  logarithmic  function 
suggested  by  Moritz  (1976,  p.  29)  appears  to  be  the  best  model  overall, 
particularly  when  the  ftmction  attains  negative  values.  The  fact  that  the  2-D 
cross-covariance  ftmction  is  not  symmetrical  complicates  the  generation  of  the 
surface  by  this  method.  It  is  possible  to  overcome  this  problem  by  using  the 
symmetrical  2-0  function  to  approximate  the  cross-covariance  surface. 

The  ideal  ftmction  would  enable  the  generation  of  all  auto  and  cross- 
covariance fimctlons  knowing  the  pertinent  parameters  of  (say)  the  anomalous 
gravity  field.  The  third-order  Markov  ftmction  suggested  by  Jordan  (1972)  has 
this  capability.  Unfortunately,  the  theoretical  relationships  did  not  agree  with 
the  actual  relationships  in  this  instance.  It  is  felt  that  this  is  an  Important  area 
for  further  research,  if  the  usefulness  of  the  2-D  covariance  ftmction  is  to  be 
fUlly  exploited. 

The  2-0  covariance  ftmction  is  capable  of  producing  results  superior 
to  those  obtained  by  the  general  function  when  certain  conditions  are  present. 

These  conditions  will  produce  large  differences  between  elements  of  the 
covariance  matrices  derived  from  the  general  covariance  analysis  and  from  the 
2-D  covariance  analysis.  They  will  occur  when: 
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(I)  anisotropic  effects  are  present  and,  because  of  the  distribution 
of  the  data,  predictions  must  be  performed  over  large  separations 
and  in  an  asymmetric  configuration,  or 

(ii)  anisotropy  is  strongly  evident  and  homogeneous  throughout  the 
field.  Such  an  effect  can  be  seen  In  areas  where  geoidal  slopes 
are  uniformly  and  consistently  large  (e.  g.  the  geoid  slope 
across  Australia).  In  feet,  under  these  conditions  the  solution 
using  die  general  function  appears  to  break  down. 

In  any  case,  a 2-D  covariance  analysis  should  be  performed  on  data 
which  shows  anisotropic  tendencies.  This  will  indicate  the  extent  of  the  azimuth 
dependence  of  the  covariance  function  and  enable  remedial  action  to  be  taken 
(e.  g.  in  the  configuration  of  the  data  used  in  subsequent  computation)  if  this 
appears  warranted. 

The  2-D  covariance  surface  may  also  provide  useful  information 
concerning  a suitable  "trend  surface"  to  be  fitted  to  the  original  data.  Knowing 
that  the  residuals  of  the  actual  data  from  the  trend  surface  should  be  isotropic, 
it  should  be  possible  to  discover  what  nature  of  surface  must  be  fitted  in  order 
to  transform  the  2-D  covariance  surface  to  a surface  of  revolution.  (This  may 
be  best  performed  in  the  spectral  domain.)  The  residuals  can  then  be  used  in  the 
stochastic  processes  with  the  knowledge  that  they  do,  in  reality  possess  isotropic 
characteristics. 
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Appendix  A1 


The  angular  frequency  (u)  and  the  frequency  in  cycles  per  second  are 
related  by: 


When  operating  in  die  spatial  domain  the  time  (? ) in  equations  (3. 2)  to 
(3.4)  are  replaced  by  the  length  r.  The  fundamental  wavelength  of  the  signal  is 
assumed  to  correspond  to  mAx  where  m,  Ax  are  defined  in  equation  (3. 5).  (In 
the  case  of  data  set  1,  it  is  seen  that  the  fimdamental  wavelength  equals  6 grid 
units,  l.e.  m-6.) 

The  frequency  associated  with  the  step  value  r (r  also  represents  the 
harmonic  of  the  fundamental  wavelength)  can  therefore  be  expressed  as: 


fr  » 


r 

2m  Ax 
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